is.inrange<-function(x,xrange){
  return((x>=min(xrange,na.rm=TRUE))&(x<=max(xrange,na.rm=TRUE)))
}

checkdistr.guessstart<-function(x,tdis){
  out<-NULL
  if(tdis=="normal"){
    out<-list(mean=mean(x,na.rm=TRUE),sd=sd(x,na.rm=TRUE))
  }
  if(tdis=="cauchy"){
    out<-list(location=as.numeric(median(x,na.rm=TRUE)),scale=as.numeric((diff(quantile(x,probs=c(0.75,0.25),na.rm=TRUE)))/2))
  }
  if(tdis=="logistic"){
    out<-list(location=as.numeric(mean(x,na.rm=TRUE)),scale=as.numeric(sqrt(var(x,na.rm=TRUE)*3/pi^2)))
  }
  if(tdis=="geometric"){
    out<-list(prob=min(c(1,1/mean(x,na.rm=TRUE)),na.rm=TRUE))
  }
  if(tdis=="exponential"){
    out<-list(rate=0.693/median(x,na.rm=TRUE))
  }
  if(tdis=="gamma"){
    trate<-mean(x,na.rm=TRUE)/var(x,na.rm=TRUE)
    out<-list(shape=mean(x,na.rm=TRUE)*trate,rate=trate)
  }
  if(tdis=="lognormal"){
    out<-list(meanlog=mean(log(x),na.rm=TRUE),sdlog=sd(log(x),na.rm=TRUE))
  }
  if(tdis=="beta"){
    tm<-mean(x,na.rm=TRUE)
    tv<-var(x,na.rm=TRUE)
    out<-list(shape1=tm*((tm*(1-tm))/tv - 1),shape2=(1-tm)*((tm*(1-tm))/tv - 1))
  }
  return(out)
}

checkdistr.single<-function(x,tdis="normal"){
  dfun<-switch(tolower(tdis),
              normal="norm",
              cauchy="cauchy",
              geometric="geom",
              exponential="exp",
              poisson="pois",
              logistic="logis",
              weibull="weibull",
              gamma="gamma",
              lognormal="lnorm",
              t="t",
              beta="beta")
  out<-data.frame(dist=tdis,fun=dfun,est_1=NA,est_2=NA,est_3=NA,name_1=NA,name_2=NA,name_3=NA,p.value=NA,d=NA,fit=FALSE,stringsAsFactors=FALSE)
  tf<-NA
  try(tf<-fitdistr(x,tdis))
  if(is.na(tf)){
    if(tdis%in%c("normal","cauchy","logistic","gamma","beta")){
      tfs<-checkdistr.guessstart(x,tdis)
      try(tf<-fitdistr(x,tdis,start=tfs))
      for(i in 1:min(c(2,length(tfs)),na.rm=TRUE)){
        out[,paste("est_",i,sep="")]<-tfs[i]
        out[,paste("name_",i,sep="")]<-names(tfs[i])
      }
    }
    if(tdis%in%c("lognormal")){
      tfs<-checkdistr.guessstart(x,tdis)
      for(i in 1:min(c(2,length(tfs)),na.rm=TRUE)){
        out[,paste("est_",i,sep="")]<-tfs[i]
        out[,paste("name_",i,sep="")]<-names(tfs[i])
      }
    }
  }
  
  if(!is.na(tf)){
    tcoeff<-complex(real=tf$estimate,imaginary=tf$sd)
    names(tcoeff)<-names(tf$estimate)
    if(tdis=="beta"){
      tkst<-ks.test(x,paste("p",dfun,sep=""),tf$estimate[1],tf$estimate[2])
    }else{
      tkst<-ks.test(x,paste("p",dfun,sep=""),tf$estimate)
    }
    out$fit<-TRUE
    out$p.value<-tkst$p.value
    out$d<-tkst$statistic
    for(i in 1:min(c(2,length(tcoeff)),na.rm=TRUE)){
      out[,paste("est_",i,sep="")]<-tcoeff[i]
      out[,paste("name_",i,sep="")]<-names(tcoeff[i])
    }
  }else{
    if(!is.na(out$est_1)){
      tfs<-checkdistr.guessstart(x,tdis)
      if(tdis=="beta"){
        tkst<-ks.test(x,paste("p",dfun,sep=""),unlist(tfs[1]),unlist(tfs[2]))
      }else{
        tkst<-ks.test(x,paste("p",dfun,sep=""),unlist(tfs))
      } 
      out$p.value<-tkst$p.value
      out$d<-tkst$statistic
    }
  }
  return(out)
}

checkdistr<-function(x=NULL){
out<-NULL
    if(!is.null(x)&require("MASS")){
    x<-x[!is.na(x)&!is.infinite(x)&!is.nan(x)]
    if(length(x)>5){
      out<-data.frame(distribution=c("Normal","Exponential","Poisson","Logistic","Geometric","Weibull","Gamma","Log-Normal","t","cauchy","Beta"),p.value=NA,
                      fun=c("norm","exp","pois","logis","geom","weibull","gamma","lnorm","t","cauchy","beta"))
      # Test auf Normal-Verteilung
      distrlist<-c("normal","cauchy","geometric","logistic","poisson")
      if(sum(x<0)==0)distrlist<-c(distrlist,"exponential")
      if(sum(x<=0)==0)distrlist<-c(distrlist,"weibull","gamma","lognormal","beta","t")
      out<-foreach(i=distrlist,.combine=rbind)%do%{
        return(checkdistr.single(x,i))
      }
      out<-out[order(out$d,decreasing=TRUE),]
      tf<-fitdistr(x,"normal")
      out[out$distribution=="Normal","p.value"]<-ks.test(x,"pnorm",tf$estimate)
      tf<-fitdistr(x,"cauchy")
      out[out$distribution=="cauchy","p.value"]<-ks.test(x,"pcauchy",tf$estimate)
      tf<-fitdistr(x,"geometric")
      out[out$distribution=="Geometric","p.value"]<-ks.test(x,"pgeom",tf$estimate)
      try({tf<-fitdistr(x,"logistic")
      out[out$distribution=="Logistic","p.value"]<-ks.test(x,"plogis",tf$estimate)})
      tf<-fitdistr(x,"Poisson")
      out[out$distribution=="Poisson","p.value"]<-ks.test(x,"ppois",tf$estimate)
      if(sum(x<0)==0){
        tf<-fitdistr(x,"exponential")
        out[out$distribution=="Exponential","p.value"]<-ks.test(x,"pexp",tf$estimate)
      }
      if(sum(x<=0)==0){
        try({tf<-fitdistr(x,"weibull")
        out[out$distribution=="Weibull","p.value"]<-ks.test(x,"pweibull",tf$estimate)})
        try({tf<-fitdistr(x,"gamma")
        out[out$distribution=="Gamma","p.value"]<-ks.test(x,"pgamma",tf$estimate)})
        tf<-fitdistr(x,"lognormal")
        out[out$distribution=="Log-Normal","p.value"]<-ks.test(x,"plnorm",tf$estimate)
        try({
          tf<-fitdistr(x,"t")
          out[out$distribution=="t","p.value"]<-ks.test(x,"pt",tf$estimate)
        })
      }
    }
    out<-out[order(out$p.value,decreasing=TRUE),]
    }
return(out)
}

check.normal<-function(indata=NULL){
if(!is.null(indata)){
  indata<-Re(indata[!is.na(indata)&!is.infinite(indata)])
  out<-data.frame(
    shapiro=shapiro.test(indata)$p.value,
    ad=ad.test(indata)$p.value,
    cvm=cvm.test(indata)$p.value,
    lillie=lillie.test(indata)$p.value,
    pearson=pearson.test(indata)$p.value
  )
}  
return(out)
}

create.chisq.table<-function(indata=NULL,factor=NULL,y=NULL){
out<-NULL
  if(!is.null(indata)){
    indata<-indata[!is.na(indata[,factor])&!is.na(indata[,y]),]
    if(is.logical(indata[,y])){
      out<-foreach(i=unique(indata[,factor]),.combine=rbind)%do%{
        return(data.frame(io=sum(indata[indata[,factor]==i,y],na.rm=TRUE),nio=sum(!indata[indata[,factor]==i,y],na.rm=TRUE),row.names=i))
      }
    }else{
      if(length(unique(indata[,y]))<(0.5*nrow(indata))){
        ystep<-unique(indata[,y])
        out<-foreach(i=unique(indata[,factor]),.combine=rbind)%do%{
          out2<-foreach(j=ystep,.combine=rbind)%do%{
            return(data.frame(value=nrow(indata[(indata[,factor]==i)&(indata[,y]==j),]),row.names=j))
          }
          out2<-t(out2)
          row.names(out2)<-i
          return(out2)
        }
      }
    }
  }
return(out)
}

cor.table<-function(indata=NULL,factors=NULL){
require(foreach)
out<-NULL
if(!is.null(indata)){
  factab<-foreach(i=1:(length(factors)-1),.combine=rbind)%do%{
    return(foreach(j=(i+1):length(factors),.combine=rbind)%do%{return(data.frame(a=as.character(factors[i]),b=as.character(factors[j]),c=NA))})
  }
  out<-foreach(i=1:nrow(factab),.combine=rbind)%dopar%{
    tft<-factab[i,]
    tdat<-indata[!is.na(indata[,as.character(tft$a)])&!is.na(indata[,as.character(tft$b)]),]
    tft$c<-cor(Re(tdat[,as.character(tft$a)]),Re(tdat[,as.character(tft$b)]))
    tft$d<-cor(1/(1e-12+Re(tdat[,as.character(tft$a)])),Re(tdat[,as.character(tft$b)]))
    return(tft)
  }
}
return(out)
}

calc.cp<-function(indata=NULL,spec=list(USL=0,OSL=1,Target=NULL),distr=NULL,p=0.05){
  out<-NULL
  if(!is.na(indata)){
    ind<-Re(indata)
    ind<-ind[!is.na(ind)]
    if(is.null(distr)){
      indistr<-checkdistr(ind)
      if(indistr[indistr$distribution=="Normal","p.value"]>.05){
        distr<-"Normal"
      }else{
        distr<-indistr[1,"distribution"]
      }
    }
    out<-data.frame(count=length(ind),USL=spec$USL,OSL=spec$OSL,mean=mean(ind,na.rm=TRUE),median=median(ind,na.rm=TRUE),distr=distr,Target=NA,cp=NA,cpk=NA,cpT=NA,cpkT=NA,cUSL=NA,cOSL=NA,ci_cpk=NA,ci_cpkT=NA)
    if(!is.null(spec$Target))out$Target=spec$Target
    out$pUSL<-sum(ind<spec$USL,na.rm=TRUE)/out$count
    out$pOSL<-sum(ind>spec$OSL,na.rm=TRUE)/out$count
    if(distr=="Normal"){#http://www.itl.nist.gov/div898/handbook/pmc/section1/pmc16.htm
      sdind<-sd(ind,na.rm=TRUE)
      cpku<-(out$mean-spec$USL)/(3*sdind)
      cpko<-(spec$OSL-out$mean)/(3*sdind)
      out$cpk<-min(c(cpku,cpko),na.rm=TRUE)
      out$cp<-sum(c(cpku,cpko),na.rm=TRUE)/2
      out$ci_cpk<-pnorm(1-p)*sqrt((1/(9*length(ind)))+out$cpk^2/(2*(length(ind)-1)))
      if(!is.null(spec$Target)){
        cpcor<-sqrt(1+((out$mean-spec$Target)/sdind)^2)
        out$cpT<-out$cp/cpcor
        out$cpkT<-out$cpk/cpcor
        out$ci_cpkT<-pnorm(1-p)*sqrt((1/(9*length(ind)))+out$cpkT^2/(2*(length(ind)-1)))
      }
      tf<-fitdistr(ind,"normal")
      pou<-pnorm(c(spec$USL,spec$OSL),tf$estimate[1],tf$estimate[2])
      out$cUSL<-pou[1]
      out$cOSL<-1-pou[2]
    }else{
      pin<-0.99865
      qou<-NA
      if(distr=="Gamma"){
        tf<-fitdistr(ind,"gamma")
        qou<-qgamma(c(pin,1-pin),tf$estimate[1],tf$estimate[2])
      }
      if(distr=="Geometric"){
        tf<-fitdistr(ind,"geometric")
        qou<-qgeom(c(pin,1-pin),tf$estimate[1])
      }
      if(distr=="Logistic"){
        tf<-fitdistr(ind,"logistic")
        qou<-qlogis(c(pin,1-pin),tf$estimate[1])
      }
      if(distr=="Poisson"){
        tf<-fitdistr(ind,"Poisson")
        qou<-qpois(c(pin,1-pin),tf$estimate[1])
      }
      if(distr=="Exponential"){
        tf<-fitdistr(ind,"exponential")
        qou<-qexp(c(pin,1-pin),tf$estimate[1])
      }
      if(distr=="Weibull"){
        tf<-fitdistr(ind,"weibull")
        qou<-qweibull(c(pin,1-pin),tf$estimate[1],tf$estimate[2])
      }
      if(distr=="Log-Normal"){
        tf<-fitdistr(ind,"lognormal")
        qou<-qlnorm(c(pin,1-pin),tf$estimate[1],tf$estimate[2])
      }
      if(distr=="t"){
        tf<-fitdistr(ind,"t")
        qou<-qt(c(pin,1-pin),tf$estimate[1],tf$estimate[2])
      }
      USL<-spec$USL
      OSL<-spec$OSL
      if(is.na(USL))USL<-0
      out$cp<-(OSL-USL)/diff(qou)
      out$cpk<-min(c(OSL-median(ind,na.rm=TRUE),median(ind,na.rm=TRUE)),na.rm=TRUE)*2/diff(qou)
      if(!is.null(spec$Target)){
        out$cpT<-(OSL-USL)/((diff(qou)/6)^2+(median(ind,na.rm=TRUE)-spec$Target)^2)^(1/6)
      }
    }
  }
  return(out)
}

analyse.cp.distr<-function(x=NULL,spec=list(USL=0,OSL=1,Target=NULL),p=0.05){
  require(foreach)
  out<-NULL
  if(!is.null(x)){
    cdist<-checkdistr(x)
    cdist<-cdist[!is.na(cdist$p.value),]
    out<-foreach(i=1:nrow(cdist),.combine=rbind)%do%{
      source("sixsigma.r")
      require(MASS)
      out2<-calc.cp(x,spec=spec,distr=cdist[i,"distribution"],p=p)
      out2<-cbind(cdist[i,],out2)
      return(out2)
    }
  }
  return(out)
}

remove.outlier<-function(x=NULL,p=0.05,distr="gamma"){
  tf<-fitdistr(x,distr)
  if(distr=="gamma")dfx<-dgamma(x,tf$estimate[1],tf$estimate[2])
  if(distr=="weibull")dfx<-dweibull(x,tf$estimate[1],tf$estimate[2])
  if(distr=="exponential")dfx<-dexp(x,tf$estimate[1],tf$estimate[2])
  if(distr=="lognormal")dfx<-dlnorm(x,tf$estimate[1],tf$estimate[2])
  if(distr=="geometric")dfx<-dgeom(x,tf$estimate[1])
  if(distr=="logistic")dfx<-dlogis(x,tf$estimate[1])
  if(distr=="poisson")dfx<-dpois(x,tf$estimate[1])
  if(distr=="t")dfx<-dt(x,tf$estimate[1],tf$estimate[2])
  out<-x[dfx>=(median(dfx)/length(x))]
  return(out)
}


nonnegres<-function(para=c(0,0,0),inval,y=NULL,negfac=1000){
  if(is.null(y))y<-rep(0,length(inval))
  ynew<-unlist(lapply(1:length(inval),function(j){return(sum(unlist(lapply(1:length(para),function(i){return(para[i]*(inval[j]^(i-1)))})),na.rm=TRUE))}))
  yres<-ynew-y
  yres[yres<0]<-yres[yres<0]*(-1)*negfac
  return(sqrt(sum(yres^2,na.rm=TRUE)))
}

determine.influence<-function(y,x,testdata){
  require(foreach)
  ux<-toupper(x)
  
  out<-foreach(ix=ux,.combine=rbind,.inorder=FALSE,.multicombine=TRUE,.packages=c("MASS"))%dopar%{
    ttd<-testdata[!is.na(testdata[,y])&!is.na(testdata[,ix]),]
    tlm<-rlm(as.formula(paste(y,ix,sep="~")),data=ttd)
    t2lm<-lm(as.formula(paste(y,ix,sep="~")),data=ttd)
    out<-data.frame(element=ix,type="v",sigma=summary(tlm)$sigma,df=summary(tlm)$df[2],stddev=summary(tlm)$stddev,cor=NA,t=summary(tlm)$coefficients[2,3],t2=summary(t2lm)$coefficients[2,4],rs=summary(t2lm)$r.squared)
    if(nrow(ttd)>5)out$cor<-cor(ttd[,y],ttd[,ix])
    ttd<-testdata[!is.na(testdata[,y])&!is.na(testdata[,ix])&is.inrange(testdata[,ix],quantile(testdata[,ix],na.rm=TRUE,probs=c(0.1,0.9))),]
    tlm<-rlm(as.formula(paste(y,ix,sep="~")),data=ttd)
    t2lm<-lm(as.formula(paste(y,ix,sep="~")),data=ttd)
    out2<-data.frame(element=ix,type="f",sigma=summary(tlm)$sigma,df=summary(tlm)$df[2],stddev=summary(tlm)$stddev,cor=NA,t=summary(tlm)$coefficients[2,3],t2=summary(t2lm)$coefficients[2,4],rs=summary(t2lm)$r.squared)
    if(nrow(ttd)>5)out2$cor<-cor(ttd[,y],ttd[,ix])
    return(rbind(out,out2))
  }
}
determine.interaction<-function(y,x,testdata,alpha=0.05){
  x<-sort(x)
  baseform<-as.formula(paste(y,paste(x,collapse="+"),sep="~"))
  basecomb<-combn(x,2)
  baseinter<-unlist(lapply(1:ncol(basecomb),function(i){return(paste(basecomb[,i],collapse=":"))}))
  testinter<-baseinter
  binter<-TRUE
  while(binter){
    testx<-c(x,testinter)
    trlm<-lm(as.formula(paste(y,paste(testx,collapse="+"),sep="~")),data=testdata)
    testinter<-testinter[testinter%in%row.names(summary(trlm)$coefficients)]
    xinter<-summary(trlm)$coefficients[testinter,4]
    if(sum(xinter>alpha)==0){
      binter<-FALSE
    }else{
      testinter<-names(xinter)[xinter<max(xinter,na.rm=TRUE)]
    }
  }
  testx<-c(x,testinter)
  binter<-TRUE
  while(binter){
    trlm<-lm(as.formula(paste(y,paste(testx,collapse="+"),sep="~")),data=testdata)
    testx<-testx[testx%in%row.names(summary(trlm)$coefficients)]
    xinter<-summary(trlm)$coefficients[testx,4]
    if(sum(xinter>alpha)==0){
      binter<-FALSE
    }else{
      testx<-names(xinter)[xinter<max(xinter,na.rm=TRUE)]
    }
  }
  return(testx)
}

outlier.regression<-function(y,x,testdata){
  dist<-unlist(lapply(1:nrow(testdata),function(i){return(summary(rlm(as.formula(paste(y,paste(x,collapse="+"),sep="~")),data=testdata[seq(1,nrow(testdata))!=i,]))$sigma)}))
}

ddist<-function(dist="normal",x,a,b){
  out<-NULL
  if(dist=="normal")out<-dnorm(x,a,b)
  if(dist=="exp")out<-dexp(x,a)
  if(dist=="weibull")out<-dweibull(x,a)
  return(out)
}

distregress_bound<-function(dist="normal",bound=list(lower=0,upper=1,aim=NULL),gran=100,alpha=0.05){
  
  if(is.null(bound$lower))bound$lower<-0
  opstart<-c(1,1)
  if(dist=="normal")opstart<-c(bound$aim,exp(1)*2/diff(range(c(bound$lower,bound$upper))))
  if(dist=="weibull"){
    bound$lower<-max(c(bound$lower,0.0001),na.rm=TRUE)
    if(is.null(bound$aim))bound$aim<-mean(c(bound$lower,bound$upper),na.rm=TRUE)
    opstart<-c(1/bound$aim,1/bound$aim^2)
  }
  if(is.null(bound$aim))bound$aim<-mean(c(bound$lower,bound$upper),na.rm=TRUE)
  top<-optim(opstart,function(i,dist,bound,alpha){
    daim<-abs(1-ddist(dist=dist,bound$aim,i[1],i[2]))
    dlow<-alpha-ddist(dist=dist,bound$lower,i[1],i[2])
    dup<-alpha-ddist(dist=dist,bound$upper,i[1],i[2])
    if(dlow<0)dlow<-dlow*2
    if(dup<0)dup<-dup*2
    return(daim+abs(dlow)+abs(dup))
  },dist=dist,bound=bound,alpha=alpha,method="L-BFGS-B",lower=0.001,upper=100)
}

reg.weibull<-function(maxvalue=0.08,faktor=10,alpha=0.05){
  tout<-optimize(function(i){return(abs((1-alpha)-pweibull(maxvalue,faktor,i)))},c(0,maxvalue))
  return(tout$minimum)
}
reg.gamma<-function(minvalue=0.5,maxvalue=1,aimvalue=NA,alpha=0.05){
  if(is.na(aimvalue))aimvalue<-(maxvalue+minvalue)/2
  tout2<-optimize(function(i){
    gshape<-1+aimvalue*i
    tout<-abs((alpha/2)-max(c(pgamma(minvalue,gshape,i),1-pgamma(maxvalue,gshape,i)),na.rm=TRUE))
    return(tout)
  },c(0,10000))
  dout<-data.frame(rate=tout2$minimum,shape=(1+aimvalue*tout2$minimum))
  return(dout)
}

check.outlier.distr<-function(refdata,testdata,p=0.05,distr="gamma"){
  tf<-fitdistr(refdata,distr)
  td<-c(refdata,testdata)
  if(distr=="gamma")dfx<-dgamma(td,tf$estimate[1],tf$estimate[2])
  if(distr=="weibull")dfx<-dweibull(td,tf$estimate[1],tf$estimate[2])
  if(distr=="exponential")dfx<-dexp(td,tf$estimate[1],tf$estimate[2])
  if(distr=="lognormal")dfx<-dlnorm(td,tf$estimate[1],tf$estimate[2])
  if(distr=="geometric")dfx<-dgeom(td,tf$estimate[1])
  if(distr=="logistic")dfx<-dlogis(td,tf$estimate[1])
  if(distr=="poisson")dfx<-dpois(td,tf$estimate[1])
  if(distr=="normal")dfx<-dnorm(td,tf$estimate[1],tf$estimate[2])
  if(distr=="t")dfx<-dt(td,tf$estimate[1],tf$estimate[2])
  out<-tail(dfx,1)<(median(dfx,na.rm=TRUE)/length(td))
  return(out)
}

